--- title: "7-4 Days alive and free of hospital by day 28" description: | This outcome was defined as the number of days alive and free of hospital from randomisation to 28 days, calculated as 28 minus the number of days spent in hospital. All patients dying within 28 days will be assigned zero free days. author: - name: James Totterdell affiliation: University of Sydney - name: Rob Mahar affiliation: University of Melbourne date: today freeze: true --- # Preamble ```{r} #| label: pkgs #| code-summary: Load packages library (tidyverse)library (labelled)library (kableExtra)library (cmdstanr)library (posterior)library (bayestestR)library (bayesplot)library (matrixStats)library (plotly)library (lubridate)library (ggdist)library (patchwork)theme_set (theme_classic (base_size = 10 , base_family = "Palatino" ) + theme (panel.grid = element_blank (),strip.background = element_blank ()))bayesplot_theme_set (theme_set (theme_classic (base_size = 10 , base_family = "Palatino" ) + theme (panel.grid = element_blank (),strip.background = element_blank ())))color_scheme_set ("red" )options (digits = 4 )``` ```{r} #| label: data #| code-summary: Prepare dataset :: load_all ()<- read_all_no_daily ()<- all_dat %>% filter_acs_itt () %>% transmute_model_cols_grp_aus_nz ()<- acs_itt_dat %>% filter (! is.na (out_dafh))# Concurrent enrolments for C2 <- acs_itt_dat %>% filter_concurrent_intermediate ()<- acs_itt_concurc2_dat %>% filter (! is.na (out_dafh))# Concurrent enrolments for C3 <- acs_itt_dat %>% filter_concurrent_std_aspirin ()<- acs_itt_concurc3_dat %>% filter (! is.na (out_dafh))# Concurrent enrolments for C4 <- acs_itt_dat %>% filter_concurrent_therapeutic ()<- acs_itt_concurc4_dat %>% filter (! is.na (out_dafh))``` ```{r} #| label: models #| code-summary: Load models <- cmdstan_model ("stan/ordinal/logistic_cumulative.stan" ) # No epoch or site <- cmdstan_model ("stan/ordinal/logistic_cumulative_epoch.stan" ) # Epoch only <- cmdstan_model ("stan/ordinal/logistic_cumulative_epoch_site.stan" ) # Full model ``` ```{r} #| label: helper-functions #| code-summary: Functions <- function (dat, format = "html" ) {<- dat %>% group_by (CAssignment = factor (CAssignment, levels = c ("C1" , "C2" , "C3" , "C4" ),labels = intervention_labels2 ()$ CAssignment[- 1 ])) %>% summarise (Patients = n (),Known = sum (! is.na (out_dafh)),Deaths = sprintf ("%i (%.0f%%)" , sum (D28_death, na.rm = TRUE ), 100 * mean (D28_death, na.rm = TRUE )),` DAFH, Median (Q1, Q3) ` = sprintf ("%.0f (%.0f, %.0f)" , median (out_dafh, na.rm = T), quantile (out_dafh, 0.25 , na.rm = TRUE ), quantile (out_dafh, 0.75 , na.rm = TRUE ))%>% bind_rows (%>% group_by (CAssignment = "Overall" ) %>% summarise (Patients = n (),Known = sum (! is.na (out_dafh)),Deaths = sprintf ("%i (%.0f%%)" , sum (D28_death, na.rm = TRUE ), 100 * mean (D28_death, na.rm = TRUE )),` DAFH, Median (Q1, Q3) ` = sprintf ("%.0f (%.0f, %.0f)" , median (out_dafh, na.rm = T), quantile (out_dafh, 0.25 , na.rm = TRUE ), quantile (out_dafh, 0.75 , na.rm = TRUE ))%>% rename (` Anticoagulation \n intervention ` = CAssignment)kable (format = format,align = "lrrrr" ,booktabs = TRUE ,linesep = "" %>% kable_styling (font_size = 9 ,latex_options = "HOLD_position" %>% row_spec (nrow (tdat), bold = T)<- function (vars = NULL ,beta_sd_var = NULL ,ctr = contr.orthonorm,p_mult = 2 ) {<- make_X_design (dat, vars, ctr)<- attr (X, "contrasts" )<- X[, - 1 ]attr (X, "contrasts" ) <- attX<- sum (grepl ("rand" , colnames (X)))<- c (rep (1 , nXtrt), beta_sd_var)<- dat$ epoch<- max (dat$ epoch)<- dat$ ctry_num<- max (region)<- dat$ site_num<- max (site)<- dat %>% :: count (ctry_num, site_num) %>% pull (ctry_num)<- ordered (dat$ out_dafh)<- as.integer (y_raw)<- dim (X)[1 ]<- dim (X)[2 ]<- length (unique (y_mod))list (N = N, K = K, J = unique_y, X = X, y = y_mod, y_raw = y_raw,M_region = M_region,M_site = M_site, site = site,M_epoch = M_epoch, epoch = epoch,region_by_site = region_by_site,p_par = p_mult * rep (1 / unique_y, unique_y),beta_sd = beta_sd) <- function (ctr = contr.orthonorm,save = FALSE <- make_primary_model_data (c ("inelgc3" , "agegte60" , "ctry" ), c (10 , 2.5 , 1 , 1 ), ctr)<- capture.output (<- ordmod$ sample (data = mdat,seed = 813508 ,chains = 8 ,parallel_chains = min (8 , parallel:: detectCores ()),iter_warmup = 1000 ,iter_sampling = 2500 ,refresh = 0 ,adapt_delta = 0.98 if (save) {$ save_object (file.path ("outputs" , "models" , "secondary" , "7-4.rds" )) <- as_draws_rvars (mfit$ draws ())# Add names <- dat %>% dplyr:: count (epoch, epoch_lab)<- dat %>% dplyr:: count (site_num, site)names (mdrws$ beta) <- colnames (mdat$ X)names (mdrws$ gamma_epoch) <- epoch_map$ epoch_labnames (mdrws$ gamma_site) <- site_map$ site# Convert to treatment log odds ratios $ Acon <- attr (mdat$ X, "contrasts" )$ randA %**% mdrws$ beta[grepl ("randA[0-9]" , names (mdrws$ beta))]$ Ccon <- attr (mdat$ X, "contrasts" )$ randC %**% mdrws$ beta[grepl ("randC[0-9]" , names (mdrws$ beta))]$ trtA <- mdrws$ Acon[- 1 ] - mdrws$ Acon[1 ]$ trtC <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]$ compare <- c ("Intermediate vs low" = exp (mdrws$ trtC[1 ]),"Low with aspirin vs low" = exp (mdrws$ trtC[2 ]),"Therapeutic vs low" = exp (mdrws$ trtC[3 ]),"Intermediate vs low with aspirin" = exp (mdrws$ trtC[1 ] - mdrws$ trtC[2 ]),"Intermediate vs therapeutic" = exp (mdrws$ trtC[1 ] - mdrws$ trtC[3 ]),"Low with aspirin vs therapeutic" = exp (mdrws$ trtC[2 ] - mdrws$ trtC[3 ])$ OR <- c (setNames (exp (mdrws$ trtC), c ("Intermediate" , "Low with aspirin" , "Therapeutic" )),"Ineligible aspirin" = exp (mdrws$ beta[which (grepl ("inelg" , colnames (mdat$ X)))]),"Age \u2265 60" = exp (mdrws$ beta[which (grepl ("age" , colnames (mdat$ X)))]),setNames (exp (mdrws$ beta[which (grepl ("ctry" , colnames (mdat$ X)))]), c ("AU/NZ" , "NP" ))return (list (dat = mdat, fit = mfit, drws = mdrws))``` # Outcome Derivation - missing, if day 28 mortality is unknown (`D28_PatientStatus` ) or if the patient was known to have been alive at day 28 but number of days in hospital (`D28_OutcomeTotalDaysHospitalised` ) is unknown- 0 if they died by day 28 (`D28_PatientStatus` )- 28 - min(28, `D28_OutcomeTotalDaysHospitalised` ) otherwise# Descriptive Analyses ## Overall `D28_OutcomeTotalDaysHospitalised` ) against the day of discharge from the index admission.```{r} #| label: fig-dis-day-vs-days-hospitalised #| code-summary: Compare D28_OutcomeTotalDaysHospitalised to days hospitalised per dischage date. #| fig-cap: | #| Number of days hospitalised against day of discharge from index admission. #| fig-height: 4 <- acs_itt_dat %>% filter (D28_death != 1 ) %>% :: count (D28_day = pmin (28 , D28_OutcomeTotalDaysHospitalised), DIS_day = pmin (28 , DIS_day))ggplot (pdat, aes (D28_day, DIS_day)) + geom_point (aes (size = n), shape = 21 ) + geom_abline () + scale_x_continuous ("Number of days hospitalised (D28_OutcomeTotalDaysHospitalised)" ,breaks = seq (0 , 30 , 2 )) + scale_y_continuous ("Day of discharge from index admission" ,breaks = seq (0 , 30 , 2 )) + labs (size = "Count" )``` ```{r} #| label: fig-dafh-dist #| code-summary: Overall distribution of DAFH #| fig-cap: | #| Distribution of days alive and free of hospital. #| fig-subcap: ACS-ITT #| fig-height: 4 <- acs_itt_nona_dat %>% :: count (dafh = ordered (out_dafh, levels = 0 : 27 ), .drop = F) %>% mutate (p = n / sum (n))<- ggplot (pdat, aes (dafh, p)) + geom_col () + labs (x = "Days alive and free of hospital" , y = "Proportion" <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )ggsave (file.path (path, "outcome-7-4-dist-overall.pdf" ), p, height = 2.5 , width = 6 )``` ```{r} #| label: tbl-7-4-summary #| code-summary: Summary of DAFH outcome by arm #| tbl-cap: Summary of days alive and free of hospital to day 28 by treatment group. save_tex_table (make_summary_table (acs_itt_dat, "latex" ), file.path ("outcomes" , "secondary" , "7-4-anticoagulation-summary" ))make_summary_table (acs_itt_dat)``` ## Age ```{r} #| label: fig-dafh-dist-age #| code-summary: DAFH by age #| fig-cap: | #| Distribution of days alive and free of hospital by age group #| fig-height: 4 <- acs_itt_nona_dat %>% :: count (agegrp = cut (AgeAtEntry, c (18 , seq (25 , 75 , 5 ), 100 ), include.lowest = T, right = F),dafh = fct_rev (ordered (out_dafh, levels = 0 : 27 )), .drop = F) %>% group_by (agegrp) %>% mutate (p = n / sum (n))<- pdat %>% group_by (agegrp) %>% summarise (n = sum (n))<- ggplot (pdat2, aes (agegrp, n)) + geom_col (colour = "grey40" , fill = "grey40" ) + geom_vline (xintercept = 60 , linetype = 2 ) + labs (y = "Number of \n participants" ,x = "Age at entry" ) + geom_vline (xintercept = 8.5 , linetype = 2 ) + theme (axis.text.x = element_text (hjust = 1 , angle = 45 ))<- ggplot (pdat, aes (agegrp, p)) + geom_col (aes (fill = dafh)) + labs (x = "Age" , y = "Cumulative \n proportion" ) + scale_fill_viridis_d ("DAFH" , option = "A" , direction = - 1 ) + guides (fill = guide_legend (reverse = TRUE )) + geom_vline (xintercept = 8.5 , linetype = 2 ) + theme (axis.text.x = element_text (hjust = 1 , angle = 45 )) + theme (legend.key.size = unit (0.25 , "lines" ))<- p1 | p2<- file.path ("outputs" , "figures" , "outcomes" , "secondary" )ggsave (file.path (path, "7-4-age.pdf" ), p, height = 2.5 , width = 6 )``` ## Country ```{r} #| label: fig-dafh-dist-country #| code-summary: DAFH by country #| fig-cap: | #| Distribution of days alive and free of hospital by country. #| fig-height: 4 <- acs_itt_nona_dat %>% :: count (dafh = ordered (out_dafh, levels = 0 : 27 ), .drop = F) %>% group_by (country) %>% mutate (p = n / sum (n)) %>% mutate (cp = cumsum (p)) %>% ungroup ()<- pdat %>% group_by (country) %>% summarise (n = sum (n))<- ggplot (pdat2, aes (country, n)) + geom_col () + labs (y = "Number of \n participants" , x = "Country of enrolment" )<- ggplot (pdat, aes (country, p)) + geom_col (aes (fill = fct_rev (dafh))) + labs (x = "Anticoagulation intervention" , y = "Cumulative \n proportion" ) + scale_fill_viridis_d ("DAFH" , option = "A" , direction = - 1 ) + guides (fill = guide_legend (reverse = TRUE )) + theme (legend.key.size = unit (0.25 , "lines" ))<- p1 | p2<- file.path ("outputs" , "figures" , "outcomes" , "secondary" )ggsave (file.path (path, "7-4-country.pdf" ), p, height = 2.5 , width = 6 )``` ## Site ```{r} #| label: fig-dafh-dist-site #| code-summary: DAFH by site #| fig-cap: | #| Distribution of days alive and free of hospital by study site #| fig-height: 4 <- all_dat %>% filter_acs_itt () %>% filter (! is.na (out_dafh)) %>% :: count (Country = factor (PT_CountryName, levels = c ("India" , "Australia" , "Nepal" , "New Zealand" ),labels = c ("India" , "Australia" , "Nepal" , "New \n Zealand" )),Site = fct_infreq (Location),dafh = ordered (out_dafh, levels = 0 : 27 )) %>% complete (dafh = ordered (0 : 27 ), nesting (Country, Site), fill = list (n = 0 )) %>% group_by (Country, Site) %>% mutate (p = n / sum (n)) %>% mutate (cp = cumsum (p)) %>% ungroup () %>% mutate (Country = droplevels (Country),Site = droplevels (Site)<- pdat %>% group_by (Country, Site) %>% summarise (n = sum (n)) %>% ungroup ()<- ggplot (pdat2, aes (Site, n)) + facet_grid ( ~ Country, scales = "free_x" , space = "free_x" ) + geom_col () + labs (y = "Number of \n participants" , x = "" ) + theme (axis.text.x = element_text (angle = 45 , hjust = 1 ),panel.border = element_rect (fill = NA ))<- ggplot (pdat, aes (Site, p)) + facet_grid ( ~ Country, scales = "free_x" , space = "free_x" ) + geom_col (aes (fill = fct_rev (dafh))) + labs (x = "Anticoagulation intervention" , y = "Cumulative \n proportion" ) + scale_fill_viridis_d ("DAFH" , option = "A" , direction = - 1 ) + guides (fill = guide_legend (reverse = TRUE , ncol = 1 )) + theme (axis.text.x = element_text (hjust = 1 , angle = 45 )) + theme (legend.key.size = unit (0.25 , "lines" ))<- p1 / p2 + plot_layout (guides = 'collect' )<- file.path ("outputs" , "figures" , "outcomes" , "secondary" )ggsave (file.path (path, "7-4-country-site.pdf" ), p, height = 4 , width = 6.25 )``` ## Calendar Time ```{r} #| label: fig-dafh-dist-calendar #| code-summary: DAFH by calendar date #| fig-cap: | #| Distribution of days alive and free of hospital by calendar time. #| fig-height: 4 <- acs_itt_nona_dat %>% :: count (yr = year (RandDate), mth = month (RandDate),dafh = ordered (out_dafh, levels = 0 : 27 ), .drop = F) %>% group_by (yr, mth) %>% mutate (p = n / sum (n)) %>% mutate (cp = cumsum (p)) %>% ungroup ()<- pdat %>% group_by (yr, mth) %>% summarise (n = sum (n)) %>% ggplot (., aes (mth, n)) + facet_grid ( ~ yr, drop = T, scales = "free_x" , space = "free" ) + geom_col () + labs (y = "Number of \n participants" , x = "Calendar date (month of year)" ) + scale_x_continuous (breaks = 1 : 12 )<- ggplot (pdat, aes (mth, p)) + facet_grid ( ~ yr, drop = T, scales = "free_x" , space = "free" ) + geom_col (aes (fill = fct_rev (dafh))) + labs (x = "Calendar date (month of year)" , y = "Cumulative \n proportion" ) + scale_fill_viridis_d ("DAFH" , option = "A" , direction = - 1 ) + guides (fill = guide_legend (reverse = TRUE )) + theme (legend.key.size = unit (0.25 , "lines" )) + scale_x_continuous (breaks = 1 : 12 )<- p1 | p2<- file.path ("outputs" , "figures" , "outcomes" , "secondary" )ggsave (file.path (path, "7-4-calendar-time.pdf" ), p, height = 2 , width = 6 )``` ## Anti-coagulation ```{r} #| label: fig-dafh-dist-anticoagulation #| code-summary: DAFH by intervention #| fig-cap: | #| Distribution of days alive and free of hospital by anti-coagulation intervention. #| fig-height: 4 <- acs_itt_nona_dat %>% :: count (CAssignment = factor (CAssignment, labels = str_replace (intervention_labels ()$ CAssignment[- 1 ], "<br>" , " \n " )),dafh = ordered (out_dafh, levels = 0 : 27 ), .drop = F) %>% group_by (CAssignment) %>% mutate (p = n / sum (n)) %>% mutate (cp = cumsum (p)) %>% ungroup ()<- ggplot (pdat, aes (CAssignment, p)) + geom_col (aes (fill = fct_rev (dafh))) + labs (x = "Anticoagulation intervention" , y = "Cumulative proportion" ) + scale_fill_viridis_d ("Days alive and \n free of hospital" , option = "A" , direction = - 1 ) + guides (fill = guide_legend (reverse = TRUE )) + theme (legend.key.size = unit (0.75 , "lines" ))<- file.path ("outputs" , "figures" , "outcomes" , "secondary" )ggsave (file.path (pth, "outcome-7-4-descriptive.pdf" ), p, height = 3 , width = 6 )``` ```{r} #| label: fig-dafh-dist-anticoagulation2 #| code-summary: DAFH by intervention #| fig-cap: | #| Cumulative distribution of days alive and free of hospital by anti-coagulation intervention. #| fig-height: 4 <- ggplot (pdat, aes (dafh, cp)) + geom_step (aes (colour = CAssignment, group = CAssignment)) + labs (x = "Days alive and free of hospital" , y = "Cumulative proportion" ) + scale_colour_viridis_d ("Days alive and \n free of hospital" , option = "A" , begin = 0 , end = 0.8 ) + guides (fill = guide_legend (reverse = TRUE ))<- file.path ("outputs" , "figures" , "outcomes" , "secondary" )ggsave (file.path (pth, "outcome-7-4-descriptive2.pdf" ), p, height = 3 , width = 6 )``` ## Sample Cumulative Logits ```{r} #| code-summary: Plot sample cumulative logits #| fig-cap: Inspect sample cumulative logits. <- acs_itt_nona_dat %>% :: count (CAssignment, out_dafh) %>% complete (CAssignment, out_dafh, fill = list (n = 0 )) %>% group_by (CAssignment) %>% mutate (p = n / sum (n))<- trt_counts %>% group_by (CAssignment) %>% mutate (clogit = logit (cumsum (p))) %>% group_by (out_dafh) %>% mutate (rel_clogit = clogit - mean (clogit)) %>% filter (out_dafh != 27 )ggplot (trt_logit, aes (CAssignment, rel_clogit)) + facet_wrap ( ~ out_dafh) + geom_point () + labs (y = "Relative (to mean) sample cumulative logit" )ggplot (trt_logit, aes (out_dafh, rel_clogit)) + facet_wrap ( ~ CAssignment) + geom_point () + labs (y = "Relative (to mean) sample cumulative logit" )``` # Modelling ## Primary Model ```{r} #| label: fit-model-primary #| code-summary: Fit primary model <- fit_primary_model (acs_itt_nona_dat, save = TRUE )``` ```{r} #| label: odds-ratio-summary-table-primary-model-save #| code-summary: Odds ratio summary table #| tbl-cap: Posterior summaries for model parameters (fixed-effects). save_tex_table (odds_ratio_summary_table_rev (res$ drws$ OR, "latex" ),"outcomes/secondary/7-4-primary-model-acs-itt-summary-table" )odds_ratio_summary_table_rev (res$ drws$ OR)``` ```{r} #| code-summary: Odds ratio summary for epoch and site #| fig-cap: Summary of epoch and site posterior odds ratios. <- plot_epoch_site_terms (exp (res$ drws$ gamma_epoch),exp (res$ drws$ gamma_site),factor (res$ dat$ region_by_site, labels = c ("India" , "Australia \n New Zealand" , "Nepal" ))<- file.path ("outputs" , "figures" , "outcomes" , "secondary" , "7-4-primary-model-epoch-site-terms.pdf" )ggsave (pth, p, width = 6 , height = 4.5 )``` ```{r} #| code-summary: Posterior contrasts #| fig-cap: Summaries on comparisons. plot_or_densities (res$ drws$ compare)``` ##### Diagnostics ```{r} $ fit$ summary (c ("alpha" , "beta" , "gamma_epoch" , "gamma_site" , "tau_epoch" , "tau_site" )) %>% print (n = Inf )$ fit$ diagnostic_summary ()``` ##### Trace plots ```{r} mcmc_trace (res$ drws["beta" ])mcmc_trace (res$ drws["alpha" ])mcmc_trace (res$ drws["gamma_site" ])mcmc_trace (res$ drws["gamma_epoch" ])``` ### Posterior Predictive ```{r} <- as.integer (levels (res$ dat$ y_raw))<- res$ drws$ y_ppc<- rfun (\(x) y_raw[x])(y_ppc)<- bind_cols (acs_itt_nona_dat, tibble (y_ppc = y_ppc_raw)) %>% mutate (CAssignment = factor (levels = names (intervention_labels_short_break ()$ CAssignment),labels = intervention_labels_short_break ()$ CAssignment))<- function (grp, d = 0 : 27 ) {%>% group_by (grp = {{grp}}) %>% summarise (y_lteq = map_dbl (d, ~ mean (out_dafh <= .x)),ypp_lteq = map (d, ~ rvar_mean (y_ppc <= .x)),y_eq = map_dbl (d, ~ mean (out_dafh == .x)),ypp_eq = map (d, ~ rvar_mean (y_ppc == .x))%>% unnest (c (ypp_lteq, ypp_eq)) %>% mutate (days = d) %>% ungroup () %>% pivot_longer (: ypp_eq, names_to = c ("response" , "event" ),names_sep = "_" ,values_to = "posterior" )<- function (dat, lab = "" , xlab = "Probability" ) {ggplot (dat %>% filter (response == "ypp" , event == "eq" ), aes (y = days)) + facet_wrap ( ~ grp, nrow = 1 , scales = "free_x" ) + stat_slabinterval (aes (xdist = posterior), fatten_point = 1 ) + geom_point (data = dat %>% filter (response == "y" , event == "eq" ), aes (y = days, x = mean (posterior)),colour = "red" ,shape = 23 ) + scale_x_continuous (breaks = c (0 , 0.3 , 0.6 ),sec.axis = sec_axis (~ . , name = lab, breaks = NULL , labels = NULL )) + scale_y_continuous ("Days" ) + theme (strip.text = element_text (size = rel (0.7 )),axis.title.x = element_text (size = rel (0.7 )),axis.text.x = element_text (size = rel (0.65 )),axis.title.y = element_text (size = rel (0.75 )),axis.title.x.bottom = element_blank ())<- grp_ppc (CAssignment)<- grp_ppc (country)<- grp_ppc (epoch)<- grp_ppc (site) %>% left_join (ppc_dat %>% dplyr:: count (site, country), by = c ("grp" = "site" ))<- plot_grp_ppc (pp_C, "Anticoagulation" , "" )<- plot_grp_ppc (pp_ctry, "Country" , "" ) <- plot_grp_ppc (pp_epoch, "Epoch" , "" )<- plot_grp_ppc (%>% filter (country == "IN" ), "Sites India" , "" )<- plot_grp_ppc (%>% filter (country == "AU" ), "Sites Australia" , "" )<- plot_grp_ppc (%>% filter (country == "NP" ), "Sites Nepal" , "" )<- plot_grp_ppc (%>% filter (country == "NP" ), "Sites New Zealand" , "" )<- (p1 | p2) / p3<- p4 / p5 / (p6 | p7)<- file.path ("outputs" , "figures" , "outcomes" , "secondary" ,"7-4-primary-model-acs-itt-ppc1.pdf" )<- file.path ("outputs" , "figures" , "outcomes" , "secondary" ,"7-4-primary-model-acs-itt-ppc2.pdf" )ggsave (pth1, g1, width = 6 , height = 5 , device = cairo_pdf)ggsave (pth2, g2, width = 6 , height = 6 , device = cairo_pdf)``` ## Sensitivity: Reduced Model (No site term) ```{r} #| code-summary: Fit reduced model without site term #| eval: false <- function (dat) {<- make_primary_model_data (dat, c ("inelgc3" , "agegte60" , "ctry" ), c (10 , 2.5 , 1 , 1 ))<- capture.output (<- ordmod_epoch$ sample (data = mdat,seed = 813508 ,chains = 8 ,parallel_chains = 8 ,iter_warmup = 1000 ,iter_sampling = 2500 ,refresh = 0 ,adapt_delta = 0.98 <- as_draws_rvars (mfit$ draws ())# Add names <- dat %>% dplyr:: count (epoch, epoch_lab)names (mdrws$ beta) <- colnames (mdat$ X)names (mdrws$ gamma_epoch) <- epoch_map$ epoch_lab# Convert to treatment log odds ratios $ Acon <- attr (mdat$ X, "contrasts" )$ randA %**% mdrws$ beta[grepl ("randA[0-9]" , names (mdrws$ beta))]$ Ccon <- attr (mdat$ X, "contrasts" )$ randC %**% mdrws$ beta[grepl ("randC[0-9]" , names (mdrws$ beta))]$ trtA <- mdrws$ Acon[- 1 ] - mdrws$ Acon[1 ]$ trtC <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]$ OR <- c (setNames (exp (mdrws$ trtC), c ("Intermediate" , "Low with aspirin" , "Therapeutic" )),"Ineligible aspirin" = exp (mdrws$ beta[which (grepl ("inelg" , colnames (mdat$ X)))]),"Age \u2265 60" = exp (mdrws$ beta[which (grepl ("age" , colnames (mdat$ X)))]),setNames (exp (mdrws$ beta[which (grepl ("ctry" , colnames (mdat$ X)))]), c ("AU/NZ" , "NP" ))return (list (dat = mdat, fit = mfit, drws = mdrws))<- fit_no_site_model (acs_itt_nona_dat)``` ```{r} #| code-summary: Odds ratio summary #| tbl-cap: Posterior summaries for model parameters (fixed-effects). #| eval: false odds_ratio_summary_table_rev (res$ drws$ OR)``` ```{r} #| code-summary: Epoch and site terms #| fig-cap: Posterior summaries for epoch terms. #| eval: false <- plot_epoch_terms (exp (res$ drws$ gamma_epoch))``` ## Sensitivity: Treatment coding ```{r} #| code-summary: Fit primary model using treatment contrasts <- fit_primary_model (acs_itt_nona_dat, ctr = contr.treatment, save = FALSE )``` ```{r} #| code-summary: Odds ratio table #| tbl-cap: Posterior summaries for model parameters. odds_ratio_summary_table_rev (res$ drws$ OR)``` ```{r} #| code-summary: Epoch and site terms #| label: odds-ratio-summary-primary-model-epoch-site-terms-prior-sens-trt-ctr #| fig-cap: | #| Posterior summaries for epoch and site terms under #| treatment coding for interventions. <- plot_epoch_site_terms (exp (res$ drws$ gamma_epoch),exp (res$ drws$ gamma_site),factor (res$ dat$ region_by_site, labels = c ("India" , "Australia \n New Zealand" , "Nepal" ))``` ## Sensitivity: Concurrent controls ### Intermediate dose - **Set**: ACS-ITT-intermediate- **Covariates**: anticoagulation intervention, age group, region of enrolment.```{r} #| code-summary: Descriptive summary table save_tex_table (make_summary_table (acs_itt_concurc2_dat, "latex" ), file.path ("outcomes" , "secondary" , "7-4-anticoagulation-concurrent-intermediate-summary" ))make_summary_table (acs_itt_concurc2_dat)``` ```{r} #| code-summary: Fit model <- function (fit = TRUE ) {<- contr.orthonorm<- model.matrix ( ~ CAssignment + agegte60 + ctry, data = acs_itt_concurc2_nona_dat, contrasts.arg = list (CAssignment = ctr))[, - 1 ]<- ordered (acs_itt_concurc2_nona_dat$ out_dafh)<- as.integer (y_raw)<- dim (X)[1 ]<- dim (X)[2 ]<- length (unique (y_mod))<- list (N = N, K = K, J = unique_y, X = X, y = y_mod, y_raw = y_raw,p_par = 2 * rep (1 / unique_y, unique_y),beta_sd = c (1 , 2.5 , 1 , 1 )if (fit) {<- capture.output (<- ordmod0$ sample (data = mdat,seed = 75136 ,chains = 8 ,parallel_chains = min (8 , parallel:: detectCores ()),iter_warmup = 1000 ,iter_sampling = 2500 ,refresh = 0 $ save_object (file.path ("outputs" , "models" , "secondary" , "7-4-dafh" , "acs-itt-intermediate.rds" ))else {<- readRDS (file.path ("outputs" , "models" , "secondary" , "7-4-dafh" , "acs-itt-intermediate.rds" ))<- as_draws_rvars (mfit$ draws ())names (mdrws$ beta) <- colnames (X)$ Ccon <- as.vector (ctr (2 ) %**% mdrws$ beta[1 ])$ Ctrt <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]$ compare <- c ("Intermediate vs low" = exp (mdrws$ Ctrt[1 ]))$ OR <- c ("Intermediate" = exp (mdrws$ Ctrt), setNames (exp (mdrws$ beta[2 : 4 ]), c ("Age \u2265 60" , "Australia/New Zealand" , "Nepal" ))) return (list (dat = mdat, fit = mfit, drws = mdrws))<- fit_acs_itt_intermediate (TRUE )<- res$ dat<- res$ fit<- res$ drws``` ```{r} #| code-summary: Sample distribution #| fig-cap: | #| Distribution of days alive and free of hospital by anti-coagulation intervention in subset. #| fig-height: 4 <- acs_itt_concurc2_nona_dat %>% :: count (CAssignment = factor (levels = c ("C1" , "C2" ),labels = c ("Low" , "Intermediate" )),dafh = ordered (out_dafh, levels = 0 : 27 ), .drop = F) %>% group_by (CAssignment) %>% mutate (p = n / sum (n)) %>% mutate (cp = cumsum (p)) %>% ungroup ()<- ggplot (pdat, aes (CAssignment, p)) + geom_col (aes (fill = fct_rev (dafh)), colour = NA ) + labs (x = "Anticoagulation intervention" , y = "Cumulative proportion" ) + scale_fill_viridis_d ("Days alive and \n free of hospital" , option = "A" , direction = - 1 ) + guides (fill = guide_legend (reverse = TRUE , ncol = 3 )) + theme (legend.key.size = unit (0.75 , "lines" ))<- file.path ("outputs" , "figures" , "outcomes" , "secondary" )ggsave (file.path (pth, "outcome-7-4-descriptive-concurrent-intermediate.pdf" ), p, height = 2.5 , width = 4 )``` ```{r} #| code-summary: Posterior contrast #| tbl-cap: Posterior densities for comparisons of interest. <- plot_or_densities (mdrws$ compare) + # geom_segment(aes(x = 1, xend = Inf, y = 1 - 0.1, yend = 1 - 0.1), # arrow = arrow(length = unit(0.5, "lines"))) + geom_text (aes (x = 1 , y = 1 - 0.1 , label = sprintf ("Pr(OR > 1) = %.2f" , Pr (RV > 1 ))), hjust = 0 , nudge_x = 0.01 ,family = "Palatino" ) + scale_x_log10 (breaks = c (0.8 , 0.9 , 1 , 1.1 , 1.2 , 1.3 , 1.4 , 1.5 , 1.6 )) + labs (y = "Comparison" )``` ```{r} #| code-summary: Odds ratio summary table #| tbl-cap: Posterior summaries for model parameters. odds_ratio_summary_table_rev (mdrws$ OR)save_tex_table (odds_ratio_summary_table_rev (mdrws$ OR, "latex" ),"outcomes/secondary/7-4-primary-model-acs-itt-concurrent-intermediate-summary-table" )``` #### Diagnostics ```{r} $ summary (c ("beta" ))$ diagnostic_summary ()``` #### Trace plots ```{r} mcmc_trace (mdrws["beta" ])mcmc_trace (mdrws["alpha" ])``` ### Low-dose with aspirin - **Set**: ACS-ITT-aspirin- **Covariates**: anticoagulation intervention, age group, region of enrolment.```{r} #| code-summary: Descriptive summary table save_tex_table (make_summary_table (acs_itt_concurc3_dat, "latex" ), file.path ("outcomes" , "secondary" , "7-4-anticoagulation-concurrent-stdaspirin-summary" ))make_summary_table (acs_itt_concurc3_dat)``` ```{r} #| code-summary: Fit model <- contr.orthonorm<- model.matrix ( ~ CAssignment + agegte60 + ctry, data = acs_itt_concurc3_nona_dat, contrasts.arg = list (CAssignment = ctr))[, - 1 ]<- ordered (acs_itt_concurc3_nona_dat$ out_dafh)<- as.integer (y_raw)<- dim (X)[1 ]<- dim (X)[2 ]<- length (unique (y_mod))<- list (N = N, K = K, J = unique_y, X = X, y = y_mod, y_raw = y_raw,p_par = 2 * rep (1 / unique_y, unique_y),beta_sd = c (1 , 1 , 2.5 , 1 )<- capture.output (<- ordmod0$ sample (data = mdat,seed = 135356 ,chains = 8 ,parallel_chains = 8 ,iter_warmup = 1000 ,iter_sampling = 2500 ,refresh = 0 <- as_draws_rvars (mfit$ draws ())names (mdrws$ beta) <- colnames (X)$ Ccon <- as.vector (ctr (3 ) %**% mdrws$ beta[1 : 2 ])$ Ctrt <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]$ compare <- c ("Intermediate vs low" = exp (mdrws$ Ctrt[1 ]),"Low with aspirin vs low" = exp (mdrws$ Ctrt[2 ]),"Intermediate vs low with aspirin" = exp (mdrws$ Ctrt[1 ] - mdrws$ Ctrt[2 ]))$ OR <- c ("Intermediate" = exp (mdrws$ Ctrt[1 ]),"Low with aspirin" = exp (mdrws$ Ctrt[2 ]),setNames (exp (mdrws$ beta[3 : 4 ]), c ("Age \u2265 60" , "Australia/New Zealand" )))``` ```{r} #| code-summary: Sample distribution #| fig-cap: | #| Distribution of days alive and free of hospital by anti-coagulation intervention in subset. #| fig-height: 4 <- acs_itt_concurc3_nona_dat %>% :: count (CAssignment = factor (levels = c ("C1" , "C2" , "C3" ),labels = c ("Low" , "Intermediate" , "Low \n with aspirin" )),dafh = ordered (out_dafh, levels = 0 : 27 ), .drop = F) %>% group_by (CAssignment) %>% mutate (p = n / sum (n)) %>% mutate (cp = cumsum (p)) %>% ungroup ()<- ggplot (pdat, aes (CAssignment, p)) + geom_col (aes (fill = fct_rev (dafh))) + labs (x = "Anticoagulation intervention" , y = "Cumulative proportion" ) + scale_fill_viridis_d ("Days alive and \n free of hospital" , option = "A" , direction = - 1 ) + guides (fill = guide_legend (reverse = TRUE , ncol = 3 )) + theme (legend.key.size = unit (0.75 , "lines" ))<- file.path ("outputs" , "figures" , "outcomes" , "secondary" )ggsave (file.path (pth, "outcome-7-4-descriptive-concurrent-stdaspirin.pdf" ), p, height = 2.5 , width = 6 )``` ```{r} #| code-summary: Posterior contrast #| tbl-cap: Decision summaries on comparisons plot_or_densities (mdrws$ compare)``` ```{r} #| code-summary: Odds ratio summary table #| tbl-cap: Posterior summaries for model parameters. odds_ratio_summary_table_rev (mdrws$ OR)save_tex_table (odds_ratio_summary_table_rev (mdrws$ OR, "latex" ),"outcomes/secondary/7-4-primary-model-acs-itt-concurrent-stdaspirin-summary-table" )``` #### Diagnostics ```{r} $ summary (c ("beta" ))$ diagnostic_summary ()``` #### Trace plots ```{r} mcmc_trace (mdrws["beta" ])mcmc_trace (mdrws["alpha" ])``` ### Therapeutic dose - **Set**: ACS-ITT-therapeutic- **Covariates**: anticoagulation intervention, age group, region of enrolment.```{r} #| code-summary: Descriptive summary table save_tex_table (make_summary_table (acs_itt_concurc4_dat, "latex" ), file.path ("outcomes" , "secondary" , "7-4-anticoagulation-concurrent-therapeutic-summary" ))make_summary_table (acs_itt_concurc4_dat)``` ```{r} #| code-summary: Fit model <- contr.orthonorm<- model.matrix ( ~ CAssignment + agegte60 + ctry, data = acs_itt_concurc4_nona_dat, contrasts.arg = list (CAssignment = ctr))[, - 1 ]<- ordered (acs_itt_concurc4_nona_dat$ out_dafh)<- as.integer (y_raw)<- dim (X)[1 ]<- dim (X)[2 ]<- length (unique (y_mod))<- list (N = N, K = K, J = unique_y, X = X, y = y_mod, y_raw = y_raw,p_par = 2 * rep (1 / unique_y, unique_y),beta_sd = c (1 , 1 , 2.5 , 1 , 1 )<- capture.output (<- ordmod0$ sample (data = mdat,seed = 49135 ,chains = 8 ,parallel_chains = 8 ,iter_warmup = 1000 ,iter_sampling = 2500 ,refresh = 0 <- as_draws_rvars (mfit$ draws ())names (mdrws$ beta) <- colnames (X)$ Ccon <- as.vector (ctr (3 ) %**% mdrws$ beta[1 : 2 ])$ Ctrt <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]$ compare <- c ("Intermediate vs low" = exp (mdrws$ Ctrt[1 ]),"Therapeutic vs low" = exp (mdrws$ Ctrt[2 ]),"Intermediate vs therapeutic" = exp (mdrws$ Ctrt[1 ] - mdrws$ Ctrt[2 ]))$ OR <- c ("Intermediate" = exp (mdrws$ Ctrt[1 ]),"Therapeutic" = exp (mdrws$ Ctrt[2 ]),setNames (exp (mdrws$ beta[3 : 5 ]), c ("Age \u2265 60" , "India" , "Australia/New Zealand" )))``` ```{r} #| code-summary: Sample distribution #| fig-cap: | #| Distribution of days alive and free of hospital by anti-coagulation intervention in subset. #| fig-height: 4 <- acs_itt_concurc4_nona_dat %>% :: count (CAssignment = factor (levels = c ("C1" , "C2" , "C4" ),labels = c ("low" , "Intermediate" , "Therapeutic" )),dafh = ordered (out_dafh, levels = 0 : 27 ), .drop = F) %>% group_by (CAssignment) %>% mutate (p = n / sum (n)) %>% mutate (cp = cumsum (p)) %>% ungroup ()<- ggplot (pdat, aes (CAssignment, p)) + geom_col (aes (fill = fct_rev (dafh))) + labs (x = "Anticoagulation intervention" , y = "Cumulative proportion" ) + scale_fill_viridis_d ("Days alive and \n free of hospital" , option = "A" , direction = - 1 ) + guides (fill = guide_legend (reverse = TRUE , ncol = 3 )) + theme (legend.key.size = unit (0.75 , "lines" ))<- file.path ("outputs" , "figures" , "outcomes" , "secondary" )ggsave (file.path (pth, "outcome-7-4-descriptive-concurrent-therapeutic.pdf" ), p, height = 2.5 , width = 6 )``` ```{r} #| code-summary: Posterior contrasts #| tbl-cap: Decision summaries on comparisons plot_or_densities (mdrws$ compare)``` ```{r} #| code-summary: Odds ratio summary table #| tbl-cap: Posterior summaries for model parameters. save_tex_table (odds_ratio_summary_table_rev (mdrws$ OR, "latex" ),"outcomes/secondary/7-4-primary-model-acs-itt-concurrent-therapeutic-summary-table" )odds_ratio_summary_table_rev (mdrws$ OR)``` #### Diagnostics ```{r} $ summary (c ("beta" ))$ diagnostic_summary ()``` #### Trace plots ```{r} mcmc_trace (mdrws["beta" ])mcmc_trace (mdrws["alpha" ])``` 
